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We present a novel technique for remote sensing of cloud droplet size distributions. 
Polarized reflectances in the scattering angle range between 135 and 165 exhibit a 
sharply defined rainbow structure, the shape of which is determined mostly by single 
scattering properties of cloud particles, and therefore, can be modeled using the Mie 
theory. Fitting the observed rainbow with such a model (computed for a parameterized 
family of particle size distributions) has been used for cloud droplet size retrievals. We 
discovered that the relationship between the rainbow structures and the corresponding 
particle size distributions is deeper than it had been commonly understood. In fact, the 
Mie theory-derived polarized reflectance as a function of reduced scattering angle (in 
the rainbow angular range) and the (monodisperse) particle radius appears to be a 
proxy to a kernel of an integral transform (similar to the sine Fourier transform on the 
positive semi-axis). This approach, called the rainbow Fourier transform (RFT), allows 
us to accurately retrieve the shape of the droplet size distribution by the application of 
the corresponding inverse transform to the observed polarized rainbow. While the basis 
functions of the proxy-transform are not exactly orthogonal in the finite angular range, 
this procedure needs to be complemented by a simple regression technique, which 
removes the retrieval artifacts. This non-parametric approach does not require any a 
priori knowledge of the droplet size distribution functional shape and is computation- 
ally fast (no look-up tables, no fitting, computations are the same as for the forward 
modeling). 

© 2012 Elsevier Ltd. All rights reserved. 


1. Introduction 


Cloud droplet size retrievals from measurements in 
the solar spectral domain can be performed using both 
total and polarized reflectances, defined respectively as 


D nl nQ_ 

R= — - and R p = -, 

/ijo H s Io 


(1) 


where I 0 is the extraterrestrial solar irradiance and /i s is 
the cosine of the solar zenith angle (SZA). The Stokes 
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parameter Q in Eq. (1) is defined with respect to the 
scattering plane containing both solar and view direc- 
tions. In the case of water cloud, the parameter U in this 
plane is 2-3 orders of magnitude smaller than Q (cf. [1]), 
and can be neglected. This allows us to use Eq. (1) as the 
definition of polarized reflectance based on the signed 
degree of linear polarization (cf. [2]). Retrievals of the 
cloud droplet size from polarized reflectance measure- 
ments in the rainbow angular range between 135° and 
165° [3,4,1] are almost free of uncertainties due to the 3D 
nature of radiation and to gaseous and aerosol absorp- 
tions, which affect the remote sensing methods based 
solely on the multispectral measurements (not including 
polarization, cf. [5,6]). This advantage is due to the fact 
that the rainbow shape (Fig. 1) is dominated by single 
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Fig. 1. Examples of polarized reflectances in the rainbow scattering 
angle range. These radiative transfer simulations were made for realistic 
cloud droplet size distributions with r e ff = 7.5pm (top) and 17.5 pm 
(bottom). In both cases the gamma distribution model was used 
with v e ff = 0.01, 0.1, and 0.2, and the cloud optical depth was assumed 
to be 5. 


scattering of light by cloud particles. To account for the 
small contributions from all other factors (multiple scat- 
tering, Rayleigh scattering, aerosol extinction, ground 
surface reflectance for thin clouds, as well as effects 
caused by rotation to the scattering plane in Eq. (1)), the 
measured polarized reflectance was fit by the expression 
of the form 

R P (&) = A ■ P'“ ie, [0; n(r )] +B ■ 9 +C, (2 ) 

where 8 is the scattering angle, n(r) is the cloud droplet 
size distribution parameterized by its effective radius r eff 
and variance v eff (the gamma distribution shape is 
assumed), while A, B, and C are the empirical fitting 
parameters (the linear term in 9 can be replaced with a 
Rayleigh-like term ~ cos 2 8 [1]). The phase matrix ele- 
ments P ( M' e> are computed using the Mie theory for a grid 
of r ef f and v e g values. Alexandrov et al. [1] demonstrated 
that this retrieval approach, while being very accurate, 
has a disadvantage associated with the necessity of 
selecting an a priori functional shape of the droplet size 
distribution. If this shape is assumed to be monomodal, 
while the retrieval algorithm is applied to the data from a 


cloud with a bimodal droplet distribution (e.g., including 
drizzle droplets), the retrieved effective size will be biased 
toward that of the dominant mode (this is a common 
feature of all least square fit retrievals). 

In this study we present an alternative retrieval 
approach based on our observation that the polarized 
reflectances in the rainbow region taken as functions of 
the (reduced) scattering angle, while parameterized by 
particle size, form a basis for an approximate integral 
transform with certain similarities to the sine Fourier 
transform on the positive semi-axis. This transformation, 
which we called the rainbow Fourier transform (RFT), 
turned out to have a simple inverse (as the ordinary 
Fourier transform). While the direct RFT acting from the 
space of droplet size distributions (functions of the 
particle radius r) into the space of rainbow polarized 
reflectances (functions of the scattering angle 9) is a 
simple integration used in direct modeling, its inverse, 
acting in the opposite direction, is a powerful retrieval 
tool for remote sensing of cloud microphysical properties. 

We will describe the difference between the exact and 
“approximate” transforms and show how to construct the 
accurate inverse RFT on examples of droplet size distribu- 
tions of various shapes (such as rectangular and bimodal 
gamma distribution). We will also present some analyti- 
cal expressions derived from approximations to the Mie 
theory, which indirectly support our method. However, 
we currently have no theoretical proof of the RFT exis- 
tence based on the Mie theory, and our method is justified 
by empirical choices and is validated by numerical 
experiments. 

In practice, the proposed retrieval method is intended 
to be used for analyses of measurements made by the 
Research Scanning Polarimeter (RSP) [7-13]. This instru- 
ment is an airborne prototype for the satellite Aerosol 
Polarimetery Sensor (APS), which was built as part 
of the NASA Glory Mission [14]. The RSP measures the I, 
Q, and U components of the Stokes vector in nine spectral 
channels with center wavelengths of 410, 470, 555, 670, 
865, 960, 1590, 1880 and 2250 nm. It is a push-broom 
sensor scanning along the aircraft track within + 60 from 
nadir (starting at forward direction) and making samples 
at 0.8° intervals. Thus, each scan consists of about 150 
instantaneous measurements. The data from the actual 
RSP scans is then aggregated into “virtual” scans, each 
consisting of all reflectances (at a variety of scattering 
angles) from a single point on the ground or at the 
cloud top. In recent years the RSP has been deployed 
onboard different aircrafts during a number of field 
campaigns, making measurements for a vast variety of 
cloud scenes. 


2. Computation of the phase matrix element P 12 for a 
polydisperse cloud 

The averaged phase matrix element P 12 (0) for a droplet 
size distribution n(r), normalized by the condition 
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is given by the expression 

Jo” <7sca(r)Pu(r,d)n(r) dr 


p 12(A) = 


J£° (Tsca(r)n(r) dr 


(4) 


Here, we use the high-resolution functions P 12 (r,0) 
and the scattering cross-section a sca (r) computed accord- 
ing to the Mie theory for the given wavelength. This 
dataset was initially developed for generating look-up 
tables for RSP/APS retrievals and has an angular resolution 
of 0.2°. It was produced by averaging the Mie theory 
results over a set of very narrow size distributions 
(triangular, of width 0.1 pm spaced every 0.05 pm up to 
the maximum radius of 100 pm). For our purposes, we 
can neglect the finite width of these distributions and 
consider the dataset as a non-averaged Mie theory output. 
In this study we use Mie datasets for three wavelengths 
representing blue, red, and IR parts of the spectrum: 
410.2, 863.5, and 2265.1 nm. The respective complex 
refractive indices of water used in these computations 
are given in Table 1. 

Eq. (4) can be written as 

rOO 

P 12 (d) = / Pi 2 (r,0)n ff (r) dr, (5) 

Jo 


where n a (r) is the “scattering” distribution 


Mr) 


cr S ca(r)n(r) 

Jo” cr sca (r)n(r) dr' 


( 6 ) 


It is normalized by the same condition (3) as n(r). The 
scattering cross-section <j sca can be expressed as 

ffsca(r) = 7tr 2 Q sra (r), (7) 


where Qs Ca (r) is the scattering efficiency. In the case 
of large and practically non-absorbing cloud 
particles Q scti = Q ex[ & 2, thus, n a (r) is close to the area 
distribution 


Mr) = 


r 2 n(r) 

Jo”r 2 n(r) dr 


( 8 ) 


normalized by 



Mr) dr = 1 . 


(9) 


Let us introduce the reduced scattering angle 
y = 6-9 0 , (10) 


where 8 0 is close to the rainbow angle 9 R [15], but 
actually is an adjustment parameter. The angle y 
varies from 0 to the angular width of the rainbow region 


Table 1 

Refractive indices and rainbow angles for the three wavelengths used in 
this study. 



Wavelength (nm) 


410.2 

863.5 

2265.1 

n r 

1.3426514 

1.3275359 

1.2815182 

rii 

1.66 x 10~ 9 

3.49 x 10~ 7 

4.17 x itr 4 

e R 

139.3° 

137.1 

129.8 

So 

137.5° 

134.5° 

123.5 


( ^ 30°). However, we can formally extend this 

range to infinity. Let us define the following functions 
of y: 

P(y) = Pu(6 0 +y) (ID 

and 

F(r,y) = P u (r,e 0 +y). (12) 

In this notation 

rOO 

P(y)= / F(r,y)n a (r) dr, (13) 

Jo 

which can be interpreted as the result of an integral 
transform in the (r,y)-space mapping the area distribution 
Mr) t0 its image n a (y) = p(y). The function F(r,y ) is the 
kernel of this transform. 


3. Integral transforms on the positive semi-axis 


An integral transform in (r,y)-space, where r> 0 
and y > 0, and its inverse have the respective general 
forms 

/* OO 

f(y)= / f(r)F(r,y)w r (r) dr, (14) 

Jo 


r OO 

f(r)= / f(y)F(r,y)w 7 (y) dy, (15) 

Jo 

where w r and w y are the corresponding weighting func- 
tions. These expressions imply the orthogonality of the 
basic functions F(r,y): 

[ F(r,,y) F(r 2 ,y)w y (y ) dy = <3( '~ 1 ~^ > (16) 

J0 Wrfrj) 

and 

f F(r,y 1 )F(r,y 2 )w r (r)dr= < ^^. (17) 

Let us consider a family of basic functions of the form 
G = G(ry), which are orthogonal with unit weight w r = 
w y = 1 : 


Jo 


G(ny) C(r 2 y) dy = <5(n-r 2 ), 


(18) 


Jo 


G(ry{)G(ry 2 ) dr = <5(y, -y 2 ). 


An example of such a family is 
G(ry) = sin ry, 


(19) 


( 20 ) 


which forms the basis of the sine Fourier transform (a similar 
example can be constructed by replacing sine with cosine in 
this expression). Other families of orthogonal basic functions 
can be created by replacing r -> r“ and y -» y^ in the argument 
of G: 


F(r,y) = r v y ll G(r 0L yP). (21) 

A straightforward substitution demonstrates that these func- 
tions obey Eqs. (14)— (17) with the weights 

w r (r) = a r a_2v ^ 1 and w r (y) = fjyP~ 2 ^ . 


( 22 ) 
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4. Analytical approximations of P 12 


The amplitude scattering matrix (cf. [16]) relates 
the components of the incident electric field (propagating 
in the positive z-direction) to that of the scattered radia- 
tion: 


f exp (-ikR+ikz) ( S\ S 4 \/£[.\ 

U/ M v S3 wUI 


(23) 


Here, R is the distance from the particle (in the far-field), 
k = 2n/X is the wave number (2 is the light wavelength), 
the indices r and / represent the components of E in any 
two mutually orthogonal directions. In the case of sphe- 
rical particles, the amplitude scattering matrix is diag- 
onal: S 3 = S 4 = 0. The Stokes vectors I={/,Q,H,V) of the 
incident and scattered radiation are related by 


C^sca 

4nR 2 


PIo, 


(24) 


where P is the 4x4 phase matrix. The element P 12 
determining polarized reflectance is expressed in terms 
of the amplitude scattering matrix elements as 

Pi2 = 7 ^(|S,i 2 -|S 2 t 2 )- (25) 

k <7 sea 


An idea of the functional dependence of P 12 on the size 
parameter 


P = kr = 


27 rr 
— 


or k = 2P 


(26) 


(r is the spherical particle radius) and on the reduced 
scattering angle in the rainbow region e = 9—9 R (here 9 R is 
the wavelength-dependent rainbow angle, see Table 1) 
can be obtained from the Airy approximation [15]. 
This analytical approximation is formally valid for 
very large size parameters p > 5000 and very small 
reduced angles c<0.5°; however, as we show below, it 
can lead to an approximate formula for P 12 which is 
accurate in a wider parameter range. The expression for 
the amplitude scattering matrix elements in the Airy 
approximation is 

S , « _2e'> ^ nu j K 7 / 6 Ai(-0.369K 2 / 3 £) (27) 

VsinS 


j= 1, 2. Here, n is the real part of the refractive index of 
water, tq = 0.0381, u 2 = 0.00786, and e"? is a phase factor 
irrelevant for absolute value computation. Substituting 
Eq. (27) into Eq. (25) yields 

P, 2 « ( ^^K 1 / 3 Ai 2 [-0.369K 2 / 3 £]. (28) 

sin c7 

The functions Ai 2 (z) do not obey orthonormality relations 
suitable for our purposes (cf. [17]), so we need to use 
further approximations. We use the fact that the Airy 
function Ai(z) is close to its asymptotic approximation 
[18] for a large negative argument 

Ai ( - Z ) " Z v ^ sin @ Z3/2 + 9 (29> 


even if the argument value is moderate. This leads to 


Ai 2 (— z) 


Z -V 2 

2n 


l + sin^z 3/2 'j 


( 30 ) 


and to another approximate formula for Pi 2 : 

0 033 

Pl2 ~ sin0£V2 [1 + sin(O3 ' C£3/2)] - (31) 

As mentioned above, the Airy approximation (28) is 
formally valid only in the very close vicinity of 9 R : 
<§1, that is a few degrees. While Fig. 2 (top) 
demonstrates precisely this, it also shows that the exact 
Mie function and the Airy approximation exhibit certain 
similarity within the entire rainbow range (135°-165°). It 
is also seen in this plot that Ai 2 (x) in Eq. (28) is very 
closely approximated by its large-argument asymptotics 
leading to Eq. (31). The exception is only for its first 
quarter-period, while some deviations at substantially 
larger arguments simply indicate the need for a more 
precise asymptotic formula [18]. Based on Eq. (31) the 
following analytical approximation for P 12 can be intro- 
duced: 

P l2 « - sin(0.358K£' 605 ) + 0.l£-0.05. (32) 

Its plot is shown in comparison with the exact curve in 
Fig. 2 (bottom). Note that we are not looking for a high 
precision approximation here, but only for an adequate 
analytical expression to be used for educated guess of the 
integral transform weights. 


5. Definition of RFT 

We can take the functional form of the first term in Eq. 
(32) as a proxy for the integral transform kernel (12): 

F(r,y) ~ y~ 3/4 sin (ry 3/z ). (33) 

Here, y replaces e, giving us the freedom to choose 9 0 
different from 9 R , and the particle radius r replaces the 
size parameter k. We omit all constant factors for clarity 
and use 1.605^3/2 (which makes a rather small phase 
change in Eq. (32)). This expression coincides with Eq. 
(21), where the function G(x)ccsinx (Eq. (20)), and 

a=l, P = 2 * v = 0, /i = - 4 . (34) 

Thus, according to Eq. (22) the corresponding weighting 
functions are 

w y = y 2 and w r =l. (35) 

This allows us to formally define the rainbow Fourier 
transform (RFT) of the area size distribution n a (r) as 

pOO 

ha(y) = / n a (r)F(r,y) dr (36) 

and its inverse 

n' a (r)= [ n a (y)F(r,y)y 2 dy, (37) 

Jo 

where the integration is performed within the rainbow 
region (y max ~ 30°), and we use the exact kernels F(r,y) 
derived from the Mie theory (Eq. (12)) rather than their 
approximations. It appears that these functions are not 
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Fig. 2. Exact (derived from the Mie theory, red) and approximate dependences of the phase matrix element P 12 on the reduced scattering angle c = 0-0 R . 
Top: comparison with the Airy approximation Eq. (28) (dark blue) and the asymptotic formula Eq. (31)) (light blue). Bottom: Comparison with the 
approximation of Eq. (32) (green). All computations were made for the size parameter p = 655, a 863 nm wavelength, and Or = 137°. (For interpretation 
of the references to color in this figure legend, the reader is referred to the web version of this article.) 


precisely orthogonal, so the loop-transform (direct fol- 
lowed by inverse) is not an identity operator. This means 
that n' a ^n a , but rather 

n a( r )^Ci • n 0 (r) + C 2 + noise, (38) 

where C, and C 2 are constants, while “noise” stands for 
artifacts to be addressed in the next section. 

The choice of the wavelength-dependent angle 9 0 is 
made by keeping n' a (r) constant, where n a (r) = 0. If 9 0 is 
too small, the bottom-line of n' a (r ) tends to be convex 
(having maximum around 50 pm), while if 9 0 is too high, 
the bottom-line is concave. The values of 9 0 determined in 
this way are summarized in Table 1 with the correspond- 
ing rainbow angles shown for comparison. The values of 
9 0 appear to be 2-6° lower than the corresponding angles 
9 r . Fig. 3 presents examples of RFT’s base functions (12) 
vs. the reduced angle y for 410, 863, and 2265 nm 
wavelengths. These plots show that at 9=9 0 (y = 0) the 
base functions are close to 0, while 9 R is in the middle of 
their first half-period (as it is in the Airy approximation 
[15]). It is also clear from these plots that the curves 
corresponding to the 410 and 863 nm wavelengths have a 


more regular structure than that at 2265 nm. This may be 
caused by stronger absorption of water at the latter 
wavelength (see Table 1). Unfortunately, this lack of 
structure impairs the retrievals significantly, hence we 
will present the results only for 410 and 863 nm 
wavelengths. 

Fig. 4 shows the results of application of the loop- 
transform to two area size distributions n a (r) of different 
functional shapes: a bimodal gamma distribution [16] and 
a rectangular distribution (which is constant over a 
selected interval, and zero otherwise). The value of the 
constant C 2 in Eq. (38) is determined from the large- 
droplet range r~ 90-100 pm, where n a (r) is assumed to 
be zero. This value is subtracted from the retrieved n' a (r), 
after which the constant C, is supposed to be determined 
from normalization condition Eq. (9). However, for our 
tests, when the initial n a (r) is known, we simply scale the 
returned distribution so it has the same maximum value 
as the initial one (in the case of the rectangular distribu- 
tion the median over the distribution top is taken instead 
of the absolute maximum). The initial distributions n a (r) 
are over-plotted (in black) for comparison. We use the 
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Fig. 3. Examples of RFT’s basis functions F(r,y ) vs. the reduced angle y 
for 10 and 100 pm droplet radii and different wavelengths: 410 nm 
(top), 863 nm (middle), and 2265 nm (bottom). 


metric [19] 

A =i f |n"(r)— n Q (r)| dr (39) 

Jo 

to quantify how close the retrieved distribution is to the 
initial one. Here n"(r ) denotes the loop-RFT result of 
Eq. (38) after the correction for the constants C, and C 2 
described above. If both distributions n"(r) and n a (r) are 


normalized to unity according to Eq. (9), this metric 
represents the fraction of the droplet ensemble that is 
“misplaced” in the retrieved n"(r) relative to the initial 
distribution n a (r). A = 0 would indicate that n"(r) = n a (r), 
while A = 1 corresponds to the case when n"(r) and n a (r ) 
do not have common support. In practice, normalization 
of n"(r) can be affected by retrieval errors, causing A to 
exceed 1. However, we still regard this metric as an 
adequate measure of the retrieval accuracy, which is 
conveniently expressed in per cent. We consider the 
accuracy to be good if A < 10%. 

6. Correction for weak orthogonality artifacts 

While the loop-transform results presented in Fig. 4 
show generally good resemblance to the initial size 
distributions, some significant artifacts still contaminate 
the retrievals causing large values of A ~ 30%. We assume 
that these artifacts are caused mainly by the imperfect 
orthogonality of the RFT basis functions F(r,y). We mod- 
eled the effects of the weak orthogonality analytically (see 
Appendix A) assuming that our basis functions can be 
represented in the form 

F(r,y) = cH(r,y)+g(r,y), (40) 

where c is a constant, H(r,y) do satisfy appropriate 
orthogonality conditions, while g(r,y) is an additional 
term depending mainly on y and only weakly on r (note 
that we do not know anything more specific about these 
components). Computations presented in Appendix A 
show that these assumptions lead to the expression for 
the loop-transform result similar to Eq. (38) 

n'(r) = c 2 n a (r) + t/(r), (41) 

where the correction function ij responsible for “noise” 
can be represented as 

'7(r) = f / u (r) + f/ s (r) + <; g (r) + C, t . (42) 

Here, r\ u {r) is a universal (independent of the size distribution 
n a (r)) con'ection function, C l( is a (distribution-dependent) 
constant, while jy s (r) and t/ (r) are distribution-dependent 
corrections, which, however, are substantially smaller in 
magnitude than rj u (r). The function ty s (r) has a rapidly 
oscillating structure and can be represented in the form 
t; s (r) = BjSj(r)+ B 0 s 0 (r), (43 ) 

where B 0 and B i are (distribution-dependent) constants, 


while 



s 0 (o= r™ 

Jo 

F(r,y)y 2 dy 

(44) 

and 



si(r)= r 

Jo 

yF(r,y)y 2 dy 

(45) 


are universal functions. The function f/ g (r) is assumed to be 
slowly varying, while its specific structure is unknown. 

Based on this model consideration and results of 
numerical tests, we can design a correction procedure 
allowing the separation and removal of the artifacts 
from n' a (r), and the restoration of the original distribution 
n a (r). We should note that, as mentioned above, the total 


M.D. Alexandrov et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 1 13 (2012) 2521-2535 


2527 




0.06 


863 nm, A = 30% 

R., = 40 |im, R 2 = 70 nm 
V 1 = 0.01, V 2 = 0.01 
W, = 0.5, W 2 = 0.5 



0.04 


40 60 

PARTICLE RADIUS, (im 



40 60 

PARTICLE RADIUS, nm 


Fig. 4. Examples of RFT loop-transforms (green) of two model area size distributions (plotted in black for comparison). Left: bimodal gamma distribution 
with equally (50% each) weighted modes having 40 and 70 pm effective radius and the same effective variance of 0.01. Right: rectangular distribution 
which is constant between particle radii of 30 and 70 pm, and zero otherwise. The plots correspond to 410 (top) and 863 nm (bottom) wavelengths. (For 
interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) 


constant contribution to ri a (r) can easily be determined 
from its value at very large r and subtracted. Thus, it is 
sufficient for us to determine all the r-dependent artifact 
functions up to a constant (even if it is distribution- 
dependent), so we will use the terms “distribution-inde- 
pendent” and “universal” neglecting such additive con- 
stant dependence. Our numerical tests confirm that the 
correction function rj(r) is mostly universal, i.e., its shape 
only weakly depends on the particular area distribution 
n a (r). Thus, as the first step of the correction procedure we 
compute > 7 (r) using a default size distribution n d (r). It is 
convenient for this purpose to use the constant size 
distribution n d (r) = 1 /r max , where r ma x = 100 pm is the 
upper limit of particle radius in our Mie theory dataset 
(n d (r > r ma x) = 0 is assumed). Thus, we define 

>U(r) = n' d (r) (46) 

i.e., the result of the loop-RFT applied to n d (r) (subtraction 
of the first term from Eq. (41) is not necessary, since it is 
constant). Plots of the resulting correction functions are 
shown in Fig. 5 (left) for 410 and 863 nm wavelengths 
(these functions are (arbitrarily) normalized to satisfy 


'7d( r max) = 0). The first step of our correction procedure 
is simply the subtraction of the pre-computed function 
< 7 d (r) from the initial loop-transform result n' a (r). The 
results of this operation (complemented by a moving- 
window smoothing) applied to the distributions from 
Fig. 4 are shown in Fig. 6 demonstrating improvement 
in retrieval accuracy from A ~ 30% to 5-6%. Analytically 
these results can be represented in the form 

KiO = KiO-naC) 

= c 2 n a (r)+A> 7 s (r) + A> 7 g (r) (47) 

(where we omit constant terms), since the universal >i u (r) 
cancels from the difference. The oscillatory artifacts in 
Fig. 6 are consistent with the shape of A t] s (r), which is a 
linear combination of s 0 (r) and s^r), both exhibiting rapid 
oscillations with a frequency ~ y max , as can be seen in 
Fig. 5 (right). Note that these oscillatory patterns have a 
different shape than the artifacts caused by the truncation 
of the integration range (ringing effects, Gibbs ripples), 
whose amplitude decreases with distance from the dis- 
tribution maximum. In order to remove these artifacts, we 
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0 20 40 60 80 100 

PARTICLE RADIUS, 


Fig. 5. Plots of the universal correction functions for 410 nm (top) and 863 nm (bottom) wavelengths. Left: t/ d (r) (green), black curves depict moving 
averages to better show the oscillating structure. Right: s 0 (r) (red) and s, (r), scaled to magnitude of s 0 (blue). (For interpretation of the references to color 
in this figure legend, the reader is referred to the web version of this article.) 


numerically compute the functions s 0 (r) and Si(r) accord- 
ing to Eqs. (44) and (45), and then perform a multivariate 
linear regression on n"(r), as the second step of our clean- 
up procedure. This regression also includes an exponen- 
tial function exp(-0.07r) that improves the retrieved 
distribution shape at small r (probably, compensating 
for the unknown f/ g (r)), and uses weighting function 
r“ 5 / 2 (both are empirical guesses). Note that while s 0 (r) 
and S](r) look almost identical (after rescaling) in Fig. 5 
(right), they have enough differences in the small particle 
size range, so both of them should be used in the 
regression. The resulting size distribution is the residue 
of this regression. The examples presented in Fig. 7 
demonstrate a successful removal of the oscillating arti- 
facts present in Fig. 6 and further accuracy improvement 
to zl~3-4%. The remaining small differences with the 
initial size distributions can be attributed to the unknown 
smooth component fj g (r), which cannot be objectively 
separated from the distribution shape. 

Note that since ;y d (r) exhibits some slope, the final 
adjustment of the lower integration limit 9 0 should be made 
based on the condition that the retrieved size distribution has 
a constant bottom-line after ij d is subtracted. 


7. Retrievals in the presence of multiple scattering 

While polarized reflectance emerging from a cloud is 
largely dominated by single scattering by cloud droplets, 
it also includes a residual contribution from multiple 
scattering and other factors (Rayleigh scattering, aerosol 
extinction, etc.). This contribution can be well approxi- 
mated by a term linear in scattering angle plus a constant 
[3,4,1], as shown in Eq. (2). The additional term 

s(y) = By+C (48) 

in n„(y) generates contamination of n' a (r) having the form 
s(r) = Bs 1 (r) + Cs 0 (r), (49) 

where s 0 (r) and s^r) are defined by Eqs. (44) and (45). The 
effect of such contamination generated by s(y) = 0.1y + 0.2 
on the retrievals from Fig. 6 is shown in Fig. 8. Obviously, 
this contamination has the same shape as r/ s (r) described 
in the previous section (Eq. (43)), and is therefore 
removed during the second step of our correction proce- 
dure leading to the same final results presented in Fig. 7. 
Thus, this term alone does not present an additional 
challenge to the retrieval method. However, the retrievals 
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Fig. 6. Same as in Fig. 4, but with the corresponding correction functions i/ d (r) [Fig. 5 (left)] subtracted. In addition to that a moving (11 points=0.5 nm) 
window smoothing was applied. The plots correspond to 410 (top) and 863 nm (bottom) wavelengths. 


from the polarized reflectances having the form of Eq. (2) 
are complicated by the presence of the unknown coeffi- 
cient A. Indeed, while h a (y) is now scaled, the default 
loop-RFT result i/ d (r) should be scaled accordingly before 
the subtraction. This suggests that we should combine the 
first and second steps of the correction method described 
above into a single multivariate regression in three 
components: >7 d (r), s 0 (r), and Sj(r) (with a possible 
additional smooth component). Generally speaking, an 
increase in the number of regression components could 
make it less stable (due to increasing possibility of trade- 
offs between these components). However, in our case the 
results of the combined regression correction procedure 
are practically indistinguishable from those of the original 
method. 

In order to test the RFT-based retrieval algorithm on 
more realistic data, we applied it to the polarized reflec- 
tances simulated by the vector doubling/adding code [20], 
which has some modifications [21] compared to the 
original version [16]. The forward modeling of these 
polarized reflectances (shown in Fig. 1) was performed 
for thick (optical depth 5) plane-parallel clouds with 
realistic droplet size distributions. It is clear from Fig. 1 
that larger droplets induce more rapid oscillations in 


polarized reflectance, while increase in v eff leads to 
smoothening of the periodic features. The results of 
inverse RFT with the regression correction procedure 
applied to these polarized reflectances are shown in 
Fig. 9. The simulations were done for gamma size dis- 
tributions with effective radii r eff = 7.5 and 17.5 pm 
(right), and effective variances v eff = 0.01 , 0.1, and 0.2. 
These numbers are parameters of number size distribu- 
tions, while the corresponding area size distributions are 
plotted. The effective radius r' eff and variance i/ ff of the 
area size distribution can be expressed in terms of their 
number size distribution counterparts as 

Kff = i ^ 2 v ff and r' ef[ = r eff O +2v eff ). (50) 

The plots show quite good agreement between the retrie- 
vals and the initial droplet size distributions used in 
forward modeling, while some artifacts remain. They 
include unphysical oscillations in the case of narrow 
(u efr = 0.01) distributions, and some shifts towards smal- 
ler sizes when r eff = 17.5 pm (0.15 pm shift at u eff = 0.01, 
0.5 pm - at Ueff = 0.1, and 0.75 pm - at i' eff = 0.2). The 
latter shift, however, may be an artifact of the RT model, 
rather than of the retrieval method. The retrieval accuracy 
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appears to be better for larger values of r eff and/or n eff , 
since the RFT artifacts are stronger in the small-radius 
range. In the worst considered case of r eff = 7.5 pm and 
v e ff = 0.01 , the error is particularly large (A = 50%), while 
for v ef f >0.1 it is substantially smaller (A < 14%). In the 
case of r eff = 17.5pm, the accuracy is even better 
(A ~ 6—10%). 

8. Practical application issues 

In this paper we focus on the introduction of a new 
retrieval algorithm and description of its mathematical 
basis. Thus, we will only briefly address the issues of 
application of RFT to real data (such as measurements by 
the RSP), leaving detailed sensitivity studies to future 
publications. Some of these sensitivities, which are com- 
mon to all rainbow-based retrieval techniques are 
described in detail in another article [1]. They include 
the effects of rotation to the principal plane, uncertainties 
in aircraft attitude (pitch and crab angles), multiple 
scattering contribution (including 3D effects), and the 
presence of aerosol layer above clouds. 

The issues specific to RFT are related to angular range and 
resolution requirements for the measurements, as well as the 
information content of the retrieved droplet distribution 


shape. Our preliminary tests on simulated and real RSP data 
showed that availability of measurements in full rainbow 
scattering angle range is essential for applicability of RFT. For 
example, reduction of the upper limit of this range from 165 
to 155° significantly impairs the retrievals. This is under- 
standable, since such reduction affects orthogonality of the 
RFT’s basis functions. Degradation of measurement resolution 
is expected to primarily affect high-frequency components of 
the angular spectrum, which correspond, as we can see from 
Fig. 1, to large particles and/or narrow size distributions (as 
well as to sharp features in the size distribution shape). A 
detailed quantitative study of these effects is yet to be done. 
The results of our preliminary tests show that the effect of a 
factor of four resolution reduction (from the model’s 0.2° to 
the RSP’s 0.8°) is practically indistinguishable for area size 
distributions having monomodal gamma shape with r eff up 
to 70 pm (it causes only a few per cent increase in A). RFT 
with 1.6 ' resolution (factor of 8 degradation) works fine for 
relatively wide distributions (u eff = 0.1 -0.2), while in the 
case of a narrow distribution (y efr = 0.01) its performance 
worsens. However, even in the latter case it still can 
adequately resolve the distribution maximum for r eff values 
up to 40 pm (but with increase in A from 7% to 35%). These 
observations suggest that RFT should not be expected to 
perform well on datasets with angular resolution worse than 


M.D. Alexandrov et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 1 13 (2012) 2521-2535 


2531 


at 

cc 

n 

a 

o 

o 


at 

at 

n 

o 

o 

o 


o 


n 

a 


LU 

at 



Fig. 8. Same as in Fig. 6, but with the effects of contamination of h 0 (y) with s(y) = 0 1 y + 0 2. 


2 °, since the droplet size distributions at the cloud top (which 
dominate the airborne polarimetric measurements) are typi- 
cally narrow. 

Unlike fitting techniques, integral transforms (including 
the RFT) cannot be forced to constrain their results to a 
certain family of functions, such as statistical distributions 
(i.e., positive normalized functions). This means that while 
RFT is expected to adequately determine the functional shape 
of a droplet size distribution, presence of noise and negative 
values may impact the integration of this distribution over 
the size range, leading to inaccurate computation of its norm 
and moments. In the case of a generic distribution shape, this 
problem needs further study. However, when the size 
distribution is expected to have gamma distribution shape, 
its effective radius and variance can be derived without 
integration. Fortunately, this is the case for the droplet size 
distributions at the cloud top, which are important for 
airborne polarimetric measurements. Their tendency to have 
gamma distribution shape was confirmed by both theoretical 
model [22,23], and airborne in situ measurements (cf. [24]). 
Our own results of application of RFT to real RSP data also 
support this assertion. While a direct functional fitting can be 
used in order to derive the distribution parameters in this 
case, we propose another, computationally simpler, approach 
based on variation of the distribution in the vicinity of its 


maximum (which is usually free from RFT-generated noise). 
This approach uses three metrics of the retrieved distribution 
shape n(r): position of the distribution maximum (mode 
radius) 

Tmax — t e jf(l — 3u e ff) (51) 

and the two ratios 

p = J- and R=^. (52) 

r max tl(T max ) 


Note that none of these metrics requires proper normal- 
ization of the retrieved distribution. Then, we derive v eff , first, 
from the relation 


lnR = 



[lnp + (l— p)] 


(53) 


for some fixed value of p and the corresponding value of R 
(we found p = 0.8 to be optimal for r > r max ). After that we 
determine r eff from Eq. (51). This method produces good 
results, when applied to the distributions from Fig. 9, espe- 
cially for v eff = 0.01, when the errors in r eff and v eff are less 
than 0.1pm and 0.01, respectively. For u eff = 0.1 and 0.2, the 
corresponding errors are larger, but they are still less than 
0.5 pm in r efr and 0.1 in v eff . 
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Fig. 9. Results of the inverse RFT (with regression correction) applied to polarized reflectances simulated using a vector radiative transfer code 
(accounting for multiple scattering contribution). The simulations were made for gamma size distributions with r eff = 7.5 (left) and 17.5 pm (right), and 
v e ff = 0.01, 0.1, and 0.2 (from top to bottom). Note that these numbers are parameters of the number size distributions, while the corresponding area size 
distributions are plotted: the retrieved (green) and the original (black). (For interpretation of the references to color in this figure legend, the reader is 
referred to the web version of this article.) 


9. Conclusions 

We described a new approach (the rainbow Fourier 
transform) to the retrieval of cloud droplet size distribu- 
tions from polarized reflectances in the rainbow scatter- 
ing angle range (135°-165°), where they are dominated 


by single scattering. This method is based on the observa- 
tion that these polarized reflectances computed using 
the Mie theory for a range of (monodisperse) particle 
radii as functions of reduced scattering angle y form a 
proxy basis of an integral transform (similar to the sine 
Fourier transform or the Bessel transform on the positive 
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semi-axis). The direct transform applied to a given droplet 
size distribution coincides with the computation of the 
corresponding polarized reflectance. The inverse trans- 
form applied to this polarized reflectance as a function of 
scattering angle yields a proxy of the original size dis- 
tribution. To obtain the distribution itself, a simple 
regression-based correction procedure is applied, which 
removes retrieval artifacts from this proxy function. Our 
analytical modeling suggested that these artifacts are 
caused by the lack of orthogonality of the RFT's base 
functions. We demonstrated good performance of the 
described technique using various sample droplet size 
distributions (bimodal gamma and rectangular) at 410 
and 863 nm wavelengths marking the boundaries of the 
visible spectral range. We also showed using RT simula- 
tions that our approach works well in realistic situations, 
when the measurements are made in the presence of 
multiple scattering. 

The main advantage of the RFT-based retrieval algo- 
rithm compared to the currently used fitting methods 
[3,4,1 ] is that this technique is non-parametric, i.e., it does 
not require any a priori knowledge of the droplet size 
distribution functional shape (including the number of 
modes). This also makes our algorithm computationally 
fast, since its computations used for inversions are essen- 
tially the same as those used for the forward modeling. 
There is no fitting using look-up tables. 

We expect that the method described has a potential 
to detect drizzle in clouds along with the cloud droplet 
distribution. However, successful use of the RFT-based 
technique imposes certain requirements on the measure- 
ments, including high angular resolution (better than 2°), 
wide angular range (close to the full rainbow range, which 
can be observed in solar principle plane measurements), 
and high measurement accuracy (low instrumental 
noise). This method will be applied to the RSP data from 
a number of field campaigns. 
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Appendix A. Effects of weak orthogonality 

Let us assume that the lack of orthogonality of the RFT 
basis functions F(r,y) is caused by an additive term g(r,y), 
i.e., these functions can be represented in the form 

F(r,y) = cH(r,y)+g(r,y), (A.l) 

where c is a constant, and H(r,y) satisfy the orthogonality 
conditions (16) and (17) with the weights from Eq. (35). 
H(r,y) can be used as the kernel of another integral 
transform 

pOO 

n H (y)= / n a (r)H(r,y) dr, (A.2) 

Jo 


which, unlike RFT, is exact, i.e. 

n a (r) = [ n H (y)H(r,y)y 2 dy. 
Jo 


(A.3) 


Here, we neglect the effect of finiteness of the integration 
limit y max , which generally may cause artificial oscilla- 
tions (ringing artifacts, Gibbs ripples), since the latter do 
not show up in numerical tests. The result of the original 
RFT (F-transform, Eq. (36)) can be expressed in terms of 
this H-transform as 


n F (y) = cn H (y)+g(y), 
where 

pO O 

g(y) = / n a (r)g(r,y) dr. 
Jo 


(A.4) 


(A.5) 


Then, the result of the inverse RFT (Eq. (37)) can be 
written as 


n' a (r)= [ [ch H (y)+g(y)]F(r,y)y 2 dy 

Jo 

= c 2 n a (r)+t](r), 
where 

i1(r) = >1 0 (r)+th(r) 

with 

'7o( r > = c f n H (y)g(r,y)y 2 dy 
Jo 

and 

t/i( r ) = / g(y)F(r,y)y 2 dy 
Jo 


(A.6) 
(A. 7) 

(A.8) 

(A.9) 


describe the non-orthogonality artifacts in the loop RFT. 
Using Eq. (A.2) and changing the integration order (we 
assume all the convergence properties necessary for this), 
Eq. (A.8) can be written as 

>7o(r) = c y 2 dy g(r,y) [ dr'n a (r')H(r , ,y) 

Jo Jo 

pOO 

= / n a (r')h(r,r')dr', (A. 10) 

Jo 


where 


g(r,y)H(r',y)y 2 dy. 


(AH) 


h(r,r') = c [‘ 

Jo 

Similarly, Eq. (A.9) can be written in the following form, 
using Eq. (A.5): 

*h(r) = c [ y 2 dy F(r,y) / d r'n a (r’)g(r',y) 

Jo Jo 

poo 

= / n a (r')f(r',r) dr', (A.12) 

Jo 


where 


/(r.r 0 = P 

Jo 


g(r,y)F(r',y)y 2 dy. 


(A.13) 


The function h(r,r ') in Eq. (A.10) can be expressed in terms 
of f(r.r’) as 

h(r,r') =f(r,r')-y(r,r'), 


(A. 14) 
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where 

y(r.r') = P m “ g(r,y)g(r J ,y)y 2 dy (A.15) 

Jo 

(obviously, y is symmetric: y(r,r') = y(r',r)). In this nota- 
tion the expression for rj(r) has the following form: 

nOO 

f](r) = / n a (r')p(r,r')dr', (A.16) 

Jo 

where 

P(t.r') =/(r,r')+/(r',r)-y(r,r') (A.17) 

is symmetric function: p(r,r' ) = p(r',r). 

Our numerical tests show that the shape of ry(r) is 
almost independent of the size distribution rt a (r) (it may, 
however, include a distribution-dependent constant, 
which can always be removed by taking into account that 
n a (r) = 0 at very large r. As we show below, complete 
distribution-independence in the above sense is achieved 
if the additive term g(r,y) does not depend on particle 
radius: g = g 0 (y). This allows us to consider the r-depen- 
dent part of the term gas a small perturbation 


g(r,y) = g 0 (y)+£ g\(r.y). (A.18) 

where e is a small parameter. In this case 

/(r,r')=/ 0 (r') + £/ 1 (r,r'), (A.19) 

where 

fo(J J ) = f go(7>F(r'.y)y 2 dy (A.20) 

Jo 

and 

/,(r,r')= [ g^(r,y)F(r',y)y 2 dy. (A.21) 

Jo 

In computation of y(r,r'), we omit terms quadratic in e 
y(r.r') =y 0 + £[y, (r) +y, (r')], ( A.22 ) 

where 

y 0 = [ go(y)v 2d "j’ < A - 23 > 

Jo 

is a constant, and 

3h(r)= [ g 0 (y)g\(r,y)y 2 dy. (A.24) 

Jo 

Thus 


P(r.r')=/o(n+/o(r') + £[/',(r,r')+/ 1 (r',r)]-y 0 -£iy 1 (r)+y 1 (r')] 

(A.25) 

and 

W) = r/ u (r) + f/ s (r) + ;; g (r) + C,,, (A.26) 

where 

'7 u ( r )=/o(r)-£yi(r) (A.27) 

is a distribution-independent (universal) function 

C„=/o-yo-£yi (A.28) 

is a constant (here the bar denotes an average with n a (r')), 
while the terms 

/•OO 

t] s (r) = e ria(r')/i(r',r) dr' (A.29) 

Jo 


and 

POO 

>]g(T ) = £ y n a( r ')f i(rj') dr' (A.30) 

are the functions of particle radius, which depend on the 
size distribution shape. Note that these terms disappear 
when £ = 0 (i.e„ when g depends on y only), and jy(r) 
becomes universal (up to a constant). To understand the 
structure of ;; s (r) and // g (r), let us represent g^r.y) by its 
Taylor’s series expansion in y 

OO 

g\(r,y)=J2su(r)y i , (A.31) 

o 

which we assume to be uniformly convergent on the 
finite interval [0,y max ] for all r. Then 

OO 

h(T,f)= ^g„( r )s i (r'), (A.32) 

0 

where 

s/(r) = [ m “ y‘F(r,y)y 2 dy. (A.33) 

Jo 

Thus 

OO 

'7s( r ) = e I^£ii s i( r ) (A.34) 

0 

and 

OO 

tl g (r) = eJ2gu(r)Si- (A.35) 

o 

The functions s,(r) show rapid oscillations with frequency 
~ y max . The plots of s 0 and S] are shown in Fig. 5 (right). 
On the other hand, g^ir) are expected to be smooth slowly 
varying functions of r. This makes it difficult (if even 
possible) to separate them from the size distribution 
features. Fortunately, they appear to be small in magni- 
tude. It appears from our numerical simulations that for 
an adequate fit of n a (r) only first two terms in the 
expansion of gi(r,y) are needed for an adequate fit of 
n'(r): 

gi(A7) SB gio(r)+gii(t)y. (A.36) 

This simplifies the expression for p s (r): 

jy s (r) « £[gioSo(r)+gi,Si(r)]. (A.37) 
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